Unambiguous determination of gravitational waveforms from binary black hole mergers 
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Gravitational radiation is properly denned only at future null infinity ( but in practice it is estimated from 
data calculated at a finite radius. We have used characteristic extraction to calculate gravitational radiation at 
J + for the inspiral and merger of two equal mass non-spinning black holes. Thus we have determined the first 
unambiguous merger waveforms for this problem. The implementation is general purpose, and can be applied 
to calculate the gravitational radiation, at J + , given data at a finite radius calculated in another computation. 
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Introductory Remarks. The computation of gravitational ra- 
diation from black hole merger events has attracted consider- 
able attention, since the pioneering work by Smarr and collab- 
orators 0-01 • With the advent of ground-based laser interfer- 
ometric gravitational wave detectors, as well as the prospect 
of the Laser Interferometer Space Antenna (LISA), interest in 
the problem has considerably increased. The measurement of 
gravitational waves will soon provide an important probe of 
strong-field nonlinear gravity, the domain of many fundamen- 
tal questions in astrophysics. The sensitivity of LISA, and of 
the upcoming advanced ground-based detectors AdLIGO and 
Ad Virgo, is so high that even an error in the waveform cal- 
culation of 0.1% (in a sense made precise later) could lead to 
an incorrect interpretation of the astrophysical properties of a 
source, or of a test of general relativity. Nowadays, there are 
several codes that can produce a stable and convergent sim- 
ulation of a black hole spacetime. However, a particular dif- 
ficulty with measuring gravitational radiation arises from the 
fact that, in general relativity, it cannot be defined locally but 
is defined only at future null infinity (J7 + ), which, physically, 
is the limit that is approached by radiation moving at the speed 
of light away from an isolated source. Since numerical evolu- 
tions are normally carried out on finite domains, there is a sys- 
tematic error caused by estimating the gravitational radiation 
from fields on a worldtube at finite radius and the uncertainty 
in how it relates to measurement at J + 0]. Even if this er- 
ror is small, the expected sensitivity of AdLIGO, AdVirgo and 
LISA implies that it is important to obtain an accurate result. 

A rigorous formalism for the global measurement of grav- 
itational energy at null infinity has been in place since the 
pioneering work of Bondi, Penrose and collaborators in the 
1960s @]; and subsequently, techniques for calculating 
gravitational radiation at J + have been developed. The idea 
which we pursue here is to combine a Cauchy or "3 + 1" nu- 
merical relativity code with a characteristic code 01. Given 
astrophysical initial data, such a method has only discretiza- 
tion error 0], and a complete mathematical specification has 
been developed 0] . There have been efforts to implement this 
method, often called Cauchv characteristic extraction (CCE), 
or characteristic extraction ITol [Till . Previous work has con- 
sidered test problems rather than that of the inspiral and 
merger of two black holes or neutron stars. Also, earlier ef- 
forts have combined the Cauchy and characteristic algorithms 



within the same code. 

Here, we describe the implementation of a CCE code as 
well as results obtained from the code for the inspiral and 
merger waveform of two equal mass, non- spinning, black 
holes. The waveforms are calculated at J + , and are thus 
the first unambiguous waveforms which have been obtained 
for this problem, in the sense of being free of gauge or finite- 
radius effects. Further, the code is general purpose, in that it is 
independent of the details of the Cauchy code, requiring only 
that it prescribes the required geometrical data on a world- 
tube. Thus its application to other astrophysical problems will 
be straightforward. 

For the specific problem of a binary black hole (BBH) 
merger, we show that the waveform obtained at J + contains 
only numerical error and is gauge-invariant. We demonstrate 
second-order convergence to zero in the amplitude and phase 
differences between two CCE runs using boundary data at dif- 
ferent radii. We compare the waveform obtained at J + with a 
finite-radius extrapolated waveform, and find that the correc- 
tions introduced by CCE are visible in the groundbased detec- 
tors AdLIGO and AdVirgo, as well as the space-based LISA 
detector. 

Cauchy Evolution. The scenario we envision is an isolated 
system (perturbed single body, or gravitationally bound bi- 
nary), in a region on which the Einstein (and possibly hydro- 
dynamic) equations must be solved. A standard procedure for 
doing this is to formulate the equations as an initial-boundary- 
value, or Cauchy, problem, in which data for the 3 -metric and 
its embedding is prescribed at a given time on a closed re- 
gion of the spacetime, E t . These are evolved according to 
the Einstein equations on the interior of the domain, and ar- 
tificial conditions on the timelike boundary, dT^. The first 
stable evolutions of a binary black hole system were carried 
out by Pretorius fl2ll . Two approaches to the evolution of 
the interior equations are in use: (a) the harmonic formu- 
lation of the Einstein equations with excised black hole in- 
teriors EH; (b) the BSSNOK (see O and references 
therein) evolution system with the black holes specified as 
moving "punctures" 

For the Cauchy evolutions used here, we have followed the 
latter approach using the formulation outlined in wh . The 
spacetime is discretised using finite differences on Cartesian 
grids and Berger-Oliger mesh refinement in the neighbour- 
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hood of the black holes lfl8ll . The wave zone is discretised by 
six overlapping coordinate patches with spherical topology. 
Interior boundary data between adjacent patches are commu- 
nicated by interpolation [0j| [13] • 

The outer boundary condition on the exterior of the domain 
is given by a linear outgoing wave condition on each of the 
evolved tensor components. Importantly, through the use of 
spherical grids in the wave zone, a sufficient resolution can 
be maintained even to a distant outer boundary, reducing the 
effect of grid reflections common in mesh-refinement codes. 
The size of the evolution domain is chosen according to the 
amount of time required, T, and the location of the outer- 
most measurement sphere, Vi . Since physical as well as con- 
straint violating modes propagate with the speed of light, an 
outer boundary located at tq^ > T + T{ ensures that mea- 
surements are causally disconnected from the influence of the 
outer boundary. 

Characteristic extraction. Our implementation of CCE 
is based on the mathematical prescription given in [9], and 
here we provide only an outline. The process is illustrated 
schematically in Fig.Q] Within a Cauchy simulation that uses 
Minkowski-like coordinates (£, x, y, z), we define a worldtube 
T by x 2 + y 2 + z 2 = r 2 , and compute the lapse a, the shift 
f3 l and the 3 -metric 7^ on T, as well as their first time and ra- 
dial derivatives. This data is then decomposed into spherical 
harmonics. The Cauchy code writes this spherical harmonic 
coefficient data to file, and later the CCE code postprocesses 
the data to reconstruct the 4-metric on the inner world tube. 
In this way, the CCE code is general purpose, as it runs inde- 
pendently of whatever Cauchy code was used to generate the 
worldtube data. 

The CCE code defines angular coordinates <p A as well as 
a time coordinate u (=£) on T, and constructs outgoing null 
geodesies with affine parameter A. It then transforms the 
Cauchy 4-metric to (u, A,</> A ) coordinates, and calculates a 
surface area radial coordinate rs, making the coordinate trans- 
formation to (u,rs,(j) A ) coordinates, in order to obtain the 
Bondi- Sachs metric data in a neighbourhood of T. This pro- 
vides the inner boundary data for the characteristic evolution, 
using a Cactus i2lll implementation of the PITT null evolution 
code with square stereographic coordinate patches ll22ll . The 
characteristic code uses coordinates based on outgoing null 
cones, and so the equations remain regular when the radial 
coordinate is compactified (by rs — » z = rs/(rs + rsr))> 
and in this way J + is included on the computational grid. 
The code computes the gravitational radiation at J + as the 
Weyl component ^4. The coordinate-independent quantity ^4 
is commonly used in numerical relativity; in appropriate coor- 
dinates, it is the second time derivative of the strain measured 
by a detector. 

Binary black hole evolution. We have carried out fully rela- 
tivistic evolutions of an equal-mass non-spinning binary black 
hole inspiral and merger. The initial data parameters for the 
closely bound black hole are determined by performing a post- 
Newtonian evolution from large separation in order to deter- 
mine the momenta for low-eccentricity (quasi-circular) trajec- 
tories 1 23]. The subsequent full nonlinear numerical relativity 
evolution proceeds for approximately 8 orbits (1350M), fol- 




FIG. 1 : Schematic of the CCE algorithm with two spatial dimensions 
suppressed. Spacelike slices, E* are evolved according to the Cauchy 
evolution scheme (horizontal lines). Geometrical data is recorded on 
a world-tube, T, which is used as interior boundary data for a char- 
acteristic evolution along u = constant null surfaces, transporting 
the data to J + . The outer boundary of the Cauchy domain, <9Et is 
chosen so that it is causally disconnected from Y over the course of 
the evolution. 



lowed by the merger and ringdown lasting another 100M. 

The evolutions have been carried out at two different grid 
resolutions in order to verify the convergence of the numerical 
scheme. The grid settings for the Cauchy code are: The cen- 
tral Cartesian grid consists of 6 levels of 2:1 mesh-refinement, 
with coarse grid spacings of h = 0.96M and h = 0.64M, 
respectively. A grid of spherical topology covers the far field, 
r G [35M, 3600M]; so that, during the time period of in- 
terest, the outer boundary is causally disconnected from any 
extraction sphere (see Fig.Q]). The radial spacing is commen- 
surate with the coarse Cartesian grid at the interface, and (for 
reasons of efficiency) is gradually scaled to h = 3.84M and 
h = 2.56M at the outer boundary for the two runs. We use a 
corresponding 7V ang = 21 and N ang = 31 points in each of 
the angular directions per patch. 

Characteristic boundary data were interpolated onto world- 
tubes located at r = 100M and r = 200M, and stored in the 
form of spherical harmonic coefficients, up to £ = 8, which 
was found to be the highest resolved mode. The resolutions 
of the characteristic evolutions are set up according to the re- 
spective resolutions of the Cauchy run. We use N r = 321 
and N r = 481 radial points, with 7V ang = 51 and N ang = 76 
angular points per angular patch. The dominant i = 2, m = 2 
mode of the gravitational waveform resulting from the numer- 
ical evolution is plotted in Fig. O to be described in more de- 
tail below. 

Invariance with respect to the world tube location is demon- 
strated in Fig [3] We have considered the differences between 
waveforms at J + resulting from two independent character- 
istic evolutions using boundary data at rr = 100M and rr = 
200M, respectively, and for two resolutions, h = 0.96M 
and h = 0.64M. The difference between the results should 
be entirely due to the discretisation error, and indeed this is 
what we find. The differences converge to zero with approx- 
imately second-order accuracy, as expected for the null evo- 
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FIG. 2: Inspiral, merger and ringdown phase of ip4(£ = 2, m = 2) as 
obtained from finite radius extrapolation (red) and at J + (blue). The 
waveforms are aligned at their peaks. There is a maximum difference 
of 1.08% in the amplitude and a dephasing of 0.019 radians between 
the two waves. These differences can introduce systematic errors to 
parameter estimation of events detected at high SNR. 
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FIG. 3: Differences A^4 in the amplitude of -04 = 2,m = 2) 
between two characteristic runs using boundary data from Rr = 
100M and Rr = 200M. The red curve shows the difference at res- 
olution h = 0.96M while the blue curve shows the difference for 
h — 0.64M, scaled so as to line up for second-order convergence. 
The expected second order convergence of our code is thus demon- 
strated. 

lution code. The figure displays the differences Aip4 in the 
amplitude of the wave mode ^±{1 = 2, m = 2) for resolu- 
tion h = 0.96M and h = 0.64M scaled for second-order 
convergence. The same order of convergence is also obtained 
for higher order modes such as ^{i = 4, m = 4) (not dis- 
played). The differences between the waveforms at J + for 
resolution h = 0.64M are of order of 0.03% in amplitude 
with a dephasing of 0.002 radians. 

The CCE waves can be used to evaluate the quality of stan- 
dard finite radius measurements, extrapolated to r —> oo; pre- 
viously, this was the most accurate option available. We do 



so by finding the Weyl component ^4 relative to a radially 
oriented null tetrad lfl7l1 (w e prefer ^4 to gauge-invariant per- 
turbative methods ll24l - [26lo . We have evaluated ^4 at six radii 
(r = 280, 300, 400,500, 600, 1000M) and extrapolated. De- 
tails are given in B19I1 . and the error is estimated as 0.03% in 
amplitude and 0.003 radians in phase. 

In Fig. [2 we compare the extrapolated waveform with 
that calculated at J + via CCE. The differences between the 
two waveforms have maximum and mean values of 1.08% 
and 0.166% in amplitude, and —0.019 and —0.004 radians 
in phase, respectively. That is, for the resolutions used, the 
numerical error in the characteristic evolution (see Fig. [3]) is 
smaller by one order of magnitude than the error between the 
extrapolated and characteristic waveforms in both amplitude 
and phase. Further, we note that the estimated error in the ex- 
trapolation is itself much smaller than the actual error between 
characteristic waveform and extrapolated waveform. This in- 
dicates that the systematic error in extrapolation has, previ- 
ously, been underestimated. The correction is towards slightly 
larger amplitudes and frequencies. CCE post-processes data 
produced by a Cauchy code, and as such there is an additional 
cost. However, it is relatively small: the Cauchy run reported 
here took ~ 336 hours, and then CCE required ~ 10 hours. 

Will the small correction to waveforms introduced by CCE 
be relevant to interpreting observational data? The answer will 
depend on the signal-to-noise ratio (SNR) of the event. At low 
SNR, whether CCE or extrapolated waveforms are used as a 
template will not affect physical interpretation. This is partic- 
ularly relevant, as numerical waveforms are being constructed 
with the intention of evaluating and parametrising detector 
templates and search algorithms l27ll . and to constrain ana- 
lytic models ll28l - [3lll . Our results indicate that extrapolations 
from a finite radius can be used to construct detector templates 
well within the accuracy standards required by matched filter- 
ing algorithms. 

However, at large SNR, the differences are significant to the 
determination of the physical parameters of a model measured 
in detector data. To demonstrate this, we follow methods de- 
scribed in l32ll33ll to determine the minimum SNR needed for 
a detected signal from a merger event to lead to different pa- 
rameter estimates depending on which waveform is used as a 
template. Table[Jdisplays the results for selected masses, indi- 
cating the maximum distance at which the difference between 
the waveforms will be relevant for the given merger event. 

The difference between the waveforms is unlikely to be rel- 
evant for LIGO, (e)LIGO and Virgo. Reasonable stellar mass 
black hole merger rates are expected only for a volume en- 
compassing sources up to a distance of at least lOOMpc. Thus, 
there may well be events detected by AdLIGO and AdVirgo 
for which the difference is important. Finally, the differences 
will certainly be relevant for LISA as they will be applicable 
to any supermassive black hole merger event throughout the 
visible universe (cH~ l is the Hubble radius). 
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Detector 


Masses 


Maximum distance 


LIGO 


50M© + 50M© 


5Mpc 


(e)LIGO 


50M© + 5OM 


8Mpc 


Virgo 


5OM + 5OM 


14Mpc 


AdLIGO 


5OM + 5OM 


197Mpc 


AdVirgo 


5OM + 5OM 


177Mpc 


LISA 


1O 7 M + 1O 7 M 


> cH' 1 



TABLE I: Maximum distance at which the difference between the 
extrapolated waveform and that at J + would be significant for a 
black hole merger event. 
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